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We present the results of our numerical analysis of a "composite" model of DNA which generalizes 
a well-known elementary torsional model of Yakushevich by allowing bases to move independently 
from the backbone. The model shares with the Yakushevich model many features and results but 
it represents an improvement from both the conceptual and the phenomenological point of view. It 
provides a more realistic description of DNA and possibly a justification for the use of models which 
consider the DNA chain as uniform. It shows that the existence of solitons is a generic feature of 
the underlying nonlinear dynamics and is to a large extent independent of the detailed modelling of 
DNA. As opposite to the Yakushevich model, where it is needed to use an unphysical value for the 
torsion in order to induce the correct velocity of sound, the model we consider supports solitonic 
solutions, qualitatively and quantitatively very similar to the Yakushevich solitons, in a fully realistic 
range of all the physical parameters characterizing the DNA. 



I. INTRODUCTION 



It is now almost thirty years that the role of sohtons in fundamental DNA functions is under serious investigation 
(see d i for a long Ust of references). In particular a lot of efforts have been made lately [1, 0, [1] to study the 
phenomenon of the so-called base-Hipping, namely the complete opening of a narrow segment of DNA, a phenomenon 
which is thought to be important for fundamental processes such as replication and transcription. 

A successful elementary model for DNA rotational modes, introduced by Yakushevich [6], makes use of a double 
chain of oscillators where every node of each chain is represented by a disc which interacts with their first neighbors 
on the same chain and with the disc in front of it on the other chain. The homogeneous version of this model succeds 
in supporting the existence of topological solitons for a wide range of the physical and geometrical parameters, also 
allowing analytical solutions for a few cases, but it was shown [6| that adding enough inhomogeneities - i.e. taking 
into account the geometrical and dynamical differences between the four possible bases for a realistic DNA segment 
- causes solitons to lose energy and finally stop their motion after a few nodes on the chain. 



I 1 
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FIG. 1: A picture (extracted from [?']) showing the phenomenon of base-flipping, namely the complete opening of a DNA's base 
pair due to the rotation of (at least) on of them by an angle of tt. It is clear from the picture that in the process of rotating 
about the DNA axis both the sugar and base planes also rotate about some other axis; to simplify as much as possible our 
model we will disregard this effect and concetrate only on the rotation about DNA's axis. 
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FIG. 2: Left. A fragment of our "composite" model of the DNA's double chain: at every node of the chain there are four 
degrees of freedom, two for each base (the i}) and two for each sugar-phosphate group (the {9„^i}). Right. Detail of a 
chain node. Every base {n,i} is allowed to rotate about the Ci atom of the corresponding sugar (see fig. [3} , represented by the 
point -Bn.i, only by an angle between [— 7r/2, 7r/2] because of the physical contraint represented by the sugar pentagon, so that 
the effective phase space of the model is (S^ x [— 7r/2, 7r/2])^"; this steric constraint is implemented in the model dynamically 
through an effective potential. See sec. Illll for an evaluation of all geometrical constants that appear in the pictures above. 



Our final goal, which will not be reached within this paper, is to determine whether or not realistic segments of 
DNA support the motion of "rotational" solitons. In the composite model proposed in by separating the degrees 
of freedom of the bases from the backbone-sugar component, we produced a model where the component supporting 
the solitons, i.e. the backbone-sugar chain, is perfectly homogeneous and the bases inhomogeneities act just as a 
perturbation of a homogenous system, leaving hope for a longer life of solitons moving on it. We leave to a future 
paper the study of the profiles of solitons in inhomogencous chains and their time evolution. 

The numerical results we present in this paper support our guess that the composite model represents an important 
improvement of the simple torsional models for DNA dynamics in the homogeneous approximation: it supports the 
existence of solitons and moreover these solitons are very close to the corresponding ones in the Yakushevich model. 
As a bonus, the increase in the geometrical (and consequentely dynamical) detail turned out to be enough to allow 
the generation of the solitons corresponding to the Yakushevich ones using for the coupling constants values which 
are compatible with the physical ones [s^l • 

As a final remark, we point out that all numerical results and estimates of the geometrical and dynamical parameters 
included in this paper are an improvement and/or an update of the corrisponding ones published in 1], which is fully 
superseded by this one. 



II. MODEL & EQUATIONS 

Our model for DNA is a homogeneous double chain of coupled double pendula which is a natural generalization 
of a well known model by Yakushevich where, at every node of each chain, the whole group base-sugar-phosphate is 
represented by a single disc centered at the chain's backbone axis: we simply split the group in two distinct discs, one 
still centered about the backbone axis and representing the sugar-phosphate group and the second rigidly rotating 
about a fixed point of the former as shown in fig. [2] (note that, in order to keep the geometry of the model as simple 
as possible, we consider the chain as plane rather than helicoidal). 

Referring to fig. [U the coordinates of the two discs centers and of the extrema of bases, whose distance we use to 
determine the intensity of the bonds between the base pairs on the same chain node, are the following: 

A„,.((-1)''V2 ,0 ,nS) 

Cn,t {{-iyh/2 + (-1)'+^ (i? COS 0n,^ + (dfcs + r) C0S(6'„,, -I- ifna)) , (-1)'+^ (-RcOS 0na + [dbs + r) C0s(6'„,, -|- ipn,^)) , uS) 

Dn4{-iyh/2 + {-iy+^ {RcosBn.t + {dbs + 2r) cos(6'„,, + (-1)'+^ (i? cos6'„,j + (4^ + 2r) cos(6l„,, + nJ) 



where h = 2R + Ar + 2dbs + d^q. 
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The dynamical evolution of our mechanical system is determined by the Lagrangian L = T — V , where T is the 
kinetic energy and V is the interaction potential. The interactions that are relevant for DNA rotational dynamics are 
five: 

1. The torsion (Vt) between next neighbor sugar-phosphate groups on the same chain, representing the torsional 
elasticity of the backbone. This force is the results of complex molecular interactions at the backbone level and 
it is to be considered an "effective" term. Following the principle of keeping the potential expressions as simple 
as possible until some good reason is found to make them more complicated and keeping in mind that 9 can 
vary on the whole circle, we use for it the "physical pendulum" periodic potential 

N-l 2 

Vt^Y.Y.^t[l~ cos(A0„,O] (II.l) 

n— 1 i—1 

where A0„ ^ = 6„^i i — dn.i- Note that relative angles between next neighbor discs never get big, so that it 
would be safe also to use its harmonic approximation Vt ~ J2n=i Kt/2[Adn.i]'^ ■ 

2. The stacking (Vs) between next neighbors bases on the same chain, representing the tt — tt bonds between 
the rings that constitute the bases. This interaction is much better understood than the previous one and in 
particular it is clear that only depends on the relative displacement between next neighbor bases, going rapidly 
to zero together with the overlapping portion of their surface, e.g. like in a Morse-like potential. However, since 
also in this case the relative angle of any two bases next to each other keeps small, we simplified its expression 
by considering it a harmonic bond depending on the "xy" distance between the centers of the bases: 

Af-l 2 ^ 

Vs = ^ ^ -^Ks d1y{Cn+l,i, Cn,i) 
n—1 i—1 

= ^ ^ -Ks{dbs + r)'^< [a(cos6'„+i,i - cos9n,i) + cos{6n+i.i + - cos(6'„,i -I- ipn,i)]'^ + 

n=l z=l 



[Q;(sin6'„+i,,; - sm0nA) + sin(6'„+i,i + ipn+i,i) - sin(6'„,i + (pn,i)] > (11.2) 



where a = R/{di,s + r). 

The pairing {Vp) between bases on opposite chains, representing the ionic bonds which keep the helices together. 
This is the best understood force among the ones we are considering and, like in the stacking case, it is known 
to go rapidly to zero a few Angstrom far from the equilibrium position; in this case though the distance 
between pairs of bases does get big when a base is flipping and therefore a harmonic approximation result rather 
unphysical. The interaction can be more realistically modelled by a Morse-like potential [8], nevertheless we 
will produce profiles also in the harmonic approximation in order to compare our results with those in [gI], so 
we will consider both cases: 

^ 2 



2 

n=l 



n=l 

7^Kp{dbs + 2rfl [/3(cos6'„,i + cos6'„,2) + cos(6'„4 + (p„,i) -I- cos(6i„,2 + Vn,2) - 2 - 2(jf + 

[/3(sin6'„,i + sin6'„,2) + sin(6'„,i -I- ipn,i) + sin(6'„,2 + Vn,2)f | (11.3) 



where (3 = R/{dbs + 2r) and Kp — 2Z?/i^ (see eq. (jlll.ip ). Note that, in order to simplify the expression of the 
elongation from the equilibrium position, we made the "contact" approximation d^q — 0, i.e. we disregarded the 
inter-bases distance in the equilibrium position; we will show numerically in sec. IIVI that this does not change 
significantly the solitons profiles. The same is knonw to hold also in the Yakushevich model. 

4. The helicoidal interactions (Vh) between nucleotides, which are mediated by water filaments (Bernal- Fowler 
filaments). This is the only ingredient of our model which is reminiscent of the helicoidal structure of DNA; 
in particular we will consider those being on opposite helices at half-pitch distance, as they are near enough 
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in three-dimensional space due to the double helical geometry. As the nucleotide move, the hydrogen bonds in 
these filaments and those connecting the filaments to the nucleotides are stretched and thus resist differential 
motions of the two connected nucleotides. We will, for the sake of simplicity and also in view of the small 
energies involved, only consider filaments forming between the sugar-phosphate groups, thus only the backbone 
angles will be involved in these interactions. Recalling that the pitch of the helix corresponds to 10 bases in the 
B-DNA equilibrium configuration we set 



N 



Vh = Kh[l - COs(0„+5,i+l 



(11.4) 



where the sum i + 1 is meant modulo 2. 

5. The sugar wall (Vsw), representing an "effective" interaction which dynamically restricts the range of the base 
angles ipn,i to some interval [(p-, f+] - ruling out this way the possibility of topologically non trivial configura- 
tions for them ~ to represent the steric constraint represented by the sugars, which prevents the corresponding 
bases from doing a complete circle about them. We set these angles to (j)± — ±7r/2 and verified numerically 
that the profiles in the composite model converge to the corresponding ones by narrowing more and more the 
interval [</?_, In order not to interfere with the dynamics close to the equilibria positions the potential must 
be as flat as possible close to the zeroes of the ipn^i and must rise rather quickly when the (pn,i approach ±7r/2. 
A natural choice, which we implemented in [l[, would be to use some high power of the tangent function but 
divergences cause problem in numerical calculations so we use in this paper some high even power of the sine 
function. 



N 



= ^ ^ I<sw sin''((/?„_i) 



(II.5) 



n— 1 i—1 



where the coupling constant Kgw must be taken big enough to prevent the bases from passing through the ±7r/2 
barrier but also not so big to interfere too much with the dynamics when the ip„^i are closer to the equilibrium 
position. 



As for the kynetic energy, a straightforward calculation shows that for the sugar-phosphate group we have 

N 2 



2 

n.i 



(II.6) 



and for the bases 



N 



Ts = 



n—1 i—1 



(^^ J + 2(1 + a COS ipn,i)^nAdnA + (1 + 2q! COS ipn,i + ' 



(II.7) 



Note that by putting the bases angles identically equal to zero the Lagrangian L ~ Tt +Ts — Vt — Vg — Vp — Vh — Vsi 
of the composite model reduces, modulo the helicoidal term, exactly to the Yakushevich homogeneous Lagrangian: 

N 2 ^ N-l 2 

L{en,^, 0, en,^, 0) = 2 ^ "^'^^^ + I] I] t^* + + "^'('^''^ + ^)'^^] " COs(A0„,O] 

n=l i—1 n—1 i—1 

N 



+ > (1 + + 2rfKp {2(1 - cos0„,i) + 2(1 - cos^„,2) - [1 - cos(0„,i - 0„,2)]} 



(II., 



Now, once an initial condition for all angles is given, the Lagrangian determines completely the evolution of the 
state as the "trajectory" q{t) — {On,i{t), (pn.i{t)) extremizing the action / = J^^^-^ L{q{t),q{t))dt. The initial states 
we are interested in are those ones which give rise to motions which do not change (much) shape, i.e. which move 
"rigidly", satisfying therefore the "discrete wave" condition qn+k,i{t) = Qn,iit — kS/v), where 6 is the distance between 
two consecutive nodes. We look at these solutions as "discrete solitons" able to move on the DNA chain. 

For a discussion about the analytical properties of the model we refer the reader to the paper [2]. Even in the 
simplest cases it is impossible to find an explicit solution to the Lagrange equations and therefore the system must 
be analyzed numerically. To this end we use the fact that in our case h = ~S/v is small enough to claim that 



Qn,i(t) 



gn(t -h) - qnjt) 
-h 



[qnit + d/v)-qn{t)] 



rA„,,[g(t)] 
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R di, 2rA 4, -Itt rfj, R R rfs.. 4,. 2tt dt, R 



FIG. 3: A node of the DNA chain with a AT pair (left) and a GC pair (right). The measures of the distances are meant as the 
projections of the actual lenghts over the axis connecting the two phosphates of the same chain node. The distance between 
the two phosphates is ft = ISA for both base pairs. 

SO that the kynetic energy becomes 

n— 1 2—1 

Ts=^Y^ ^ {[A(^„,,]2 + 2(1 -^aCOSipn,^)^Vn,^^0n,^ + (1 + 2a COS + ) [A6l„,,] ^ } 

n— 1 i—1 

Having discretized time we are left with just one degree of freedom, the discrete space coordinate n, so that what was 
previously the system lagrangian L = X]^=i Si=i ^ri,i has become now the action (summed over n) of the discrete 
Lagrangian _L„ = Ln,i] the discrete Lagrange equations, by the minimum action principle, are clearly given 

simply by dL — ) = and its solutions, subjected to opportune initial conditions, are exactly the 

profiles of the solitons that we want to make sure to be present in this generalization of the Yakushevich model. 

III. PHYSICAL VALUES OF PARAMETERS & DISPERSION RELATIONS 

The evaluation of the main geometrical and dynamical physical quantities entering in our model is far from being a 
trivial matter. There are still very few direct measurements, for example, of the stacking energy between bases within 
DNA and the torsional backbone force - due to quite complex interactions between the sugar and the phosphate 
atoms - is just an effective term for which it is hard in principle to define a strenght. Also for geometrical parameters 
things are not better and often in the literature different values for the same physical quantities are used. In this 
section we are going to give an estimate for the values of the parameters of our model. 

A. Kinematical parameters 

The kinematical parameters include the geometrical parameters, the mass m and the momenta of inertia /. We 
opted to directly estimate them starting from PDB data [§| (note, as already mentioned, that the present results 
supersede the ones we used in The masses can be easily calculated on the basis of the chemical structure of 
DNA. Estimate of momenta of inertia is more involved, especially for the sugar-phosphate group, since the evaluation 
depends on the choice of the rotation axis and, as pointed out earlier, base-flipping rotations are a complex subject. In 
our evaluations we made the simplest approximation, namely momenta of inertia of sugar rings have been calculated 
with respect to rotations about the main axis of DNA and passing by P atom of the backbone (see fig. [3]) and inertia 
momenta of bases have been calculated with respect to rotations about the same axis by Ci atom of the sugar ring. 

The geometrical parameters of interest for our model are the longitudinal width of bases lb and of the sugar Ig^ 
the distances of the bases from the relative sugars dbs and the distance of a base from the relative dual base d^q- 
We give our estimates for the masses, moments of inertia and the geometrical parameters for the different bases and 
their mean values in table [H From those data and using the notations R — Ig, r — Idj^, h — 1^ + dbs + h + deq/2, 
he ^ ls + dbs + h (the bar denotes the mean value) , we get the average values for the geometrical parameters appearing 
in our Lagrangian shown in table |TT1 
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A 


T 


G 


C 


mean 


Sugar 


m 


134 


125 


150 


110 


130 


85 


I 


3.6 X 10^ 


3.0 X 10^ 


4.4 X 10^ 


2.3 X 10^ 


3.3 X 10^ 


2.9 X 10^ 


I 


3.9 


2.9 


4.1 


2.7 


3.4 


3.1 


dbs 


1.0 


1.0 


1.0 


1.0 


1.0 






3.0 


3.0 


3.0 


3.0 


3.0 





TABLE I: Order of magnitude for the basic geometrical parameters of the DNA. Units of measure are: atomic unit for masses 
m, 1.67x 10~*'^Kg-m^ for the inertia momenta /, Angstrom for I (the longitudinal width of bases and sugar), dbs (the distances 
sugar-base) and deq (the distance at the equilibrium for the pairs AT and GC). These values have been extracted from the 
sample "generic" B-DNA PDB data [3, provided by the Glactone Project [TO], and double checked with the data from pT| . 
that agrees within 5%. As shown in fig. [3] the lengths of bases and sugar were taken by projecting them on the direction 
passing through the two phosphate atoms of the chain node to best fit the geometry of the model; the (real) measure for the 
diameters of the bases A, T, G, C, and the sugar and for dts are respectively (in A): 4.6, 4.0, 5.7, 4.0, 3.3 and 1.5.) 



R 


r 


dbs 


deq 


h 




a 


P 


3.1 


1.7 


1.0 


3.0 


18 


15 


1.1 


0.7 



TABLE IL Numerical values (in A) of the geometrical parameters characterizing our model. Here h — 2R + 4r + 2dbs + deq is 
the width of DNA and he — h\dg=o is its value in the contact approximation. 

B. Dynamical parameters 

The determination of the four coupling constants {Kp, Kt, Ks and Kh) characterizing our model requires some 
assessments. Their order of magnitude can be estimated by considering the typical energies of hydrogen bonds (Kp) 
and the experimental results for the torsional rigidity of the chain (Kt and Kg). As we will show later, the geometry of 
our composite model allows - with values for coupling constants within the experimental ranges - to make predictions 
about dynamical properties of DNA, as the induced optical frequencies and the phonon velocities, of the same order 
of magnitude of those experimetally observed (conversely, dispersion relations allow to restrict the selection of the 
coupling constants to those values which fit better our model). 

1 . Pairing 

The coupling constant Kp of the pairing potential pi.3[) can be determined by considering the typical energy of 
hydrogen bonds. The pairing interaction involves two (in the A — T case) or three (in the G — C case) ionic hydrogen 
bonds. We can model the pairing potential with a Morse function like 

Vp{d) = D[l - e-A'('i-'^=,)]2 = ^Kp{d - d,,f + 0(d3) , (IILi) 

where D is the potential depth, d the distance between the points I3„,i and Dn 2 (see fig.[2|), deq the bond equilibrium 
distance and ^ a parameter that controls the width of the well. Different estimates of the parameters appearing in 
the potential pil.l|l are present in the literature. The estimates 

Dat = 0.030 eV, Dgc = 0.045 eV, ^iAT = 1.9 A"\ ^igc = 2.5 A"^ 




1 2 3 4 5 6 



FIG. 4: Comparison of Morse potential and its correspondig harmonic approximation 
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are given in [T3| and used in [T^ . The values 

D = 0.040 eV, = 4.45 
are given in and used in [3, [l^, [l^ . Finally, the estimates 

Dat = 0.050 eV, Doc = 0.075 eV, fiAx = ^^gc = 4A"^ 

are given in and used in [13, UB] ■ The values of coupling constants corresponding to these different values for the 
parameters appearing in the Morse potential range across an order of magnitude: 

3.5 N/m < Kp := 2D^^ < 38N/m . (III.2) 

In our numerical investigations we will use a value of Kp near to the lower bound given in (jIII.2p : that is, we adopt 
the value Kp = 4 N/m - corresponding to D = 0.030 eV, /i = 2 A ^ - which leads to an optical excitation threshold 
of 94 — yJIKpjmh — 32 cm~^ (see eq. (jIII.12p ). so to be in agreement with [ioj . 



2. Stacking 



The determination of the torsion and stacking coupling constants is more challenging and is based on a smaller 
amount of experimental data. The main information is the total torsional ri gidi ty of the DNA chain C = SS, where 
S — 3.4A is the base-pair spacing and S is the torsional rigidity. It is known [20|, |2l| that 

10"^^ J -m < C <4-10"^^J-m. (IIL3) 

This information is used e.g. in [l^, [1^, whose estimate is based on the evaluation of the free energy of superhelical 
winding; this fixes the range for the total torsional energy to be 

180kJ/mol < 5 < 720kJ/mol. (IIL4) 

In our composite model the total torsional energy of the DNA chain has to be considered as the sum of two parts, 
the base stacking energy and the torsional energy of the sugar-phosphate backbone. In order to extract the stacking 
coupling constant we use the further information that tt — it stacking bonds amount at the most to 50kJ / mol [13, [1^ ■ 
The stacking potential can be described by a Morse potential whose width a is of the order lA - so that (see eq. (|IL2[) ) 
the harmonic approximation is justified ~ and assuming for the stacking energy the highest value possible, the coupling 
constant value is Kg = 2Es/a'^ ~ 16.6N/m. The phonon speed induced by this value of Kg is C3 = 6y/ Kg /nib — 3km/s 
(see eq. (|III.12p ). which is rather close to the the estimate of 1.8km/s < Ci < 3.5km/s given in @. 



3. Torsion and helicoidal couplings 



Once picked up the stacking component, our estimate for the torsional coupling constant Kt falls in the range 

130kJ/mol < Kt < 670kJ/mol. (III.5) 

Assuming (see [26^) that Kh ~ Kt/25 = 5kJ/mol, so that C4 = y/2Kt/Is (see equations piI.12|) and piI.13p ). aU of 
these values for Kt induce phonon speeds slightly higher with respect to the estimates cited above, between 5 km/s and 
11 km/s. For our numerical investigations, to keep the phonon speed as low as possible, we will set Kt — 130kJ/mol. 



4- Sugar wall 



In order to have a barrier flat enough close to the equilibrium position we chose to use the exponent k = 100 in III. 51 
Numerical tests show that a value of at least Kgw = lO^kJ/mol is needed in order to keep the bases angles within 
the [— 7r/2,7r/2] interval. 

For numerical calculations it is more convenient to renormalize the Lagrangian with respect to some energy value 
in order to work with dimensionless quantities. To be coherent with the Yakushevich work we renormalized the 
Lagrangian with respect to the pairing coefficient Ip = Kp{dbs + 2r)^/2 ~ 223kJ/mol; the values of renormalized 
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Kt 










130 kJ/mol 


16.6 N/m 


3.5 N/m 


5 kJ/mol 


2 • W kJ/mol 


gt 


9s 


gp 


3h 


Qsw 


0.58 


1.62 


1 


0.02 


100 



TABLE III: Values of the coupling constants (above) and of their normalized counterpart (below) . 



coefficients for coupling constants gt = Kt/lp, gs = Ks{dbs + 1")'^ /'^/Ip, gn — K^/lp and — Ksw/lp are summarized 
in table HTH 

Comparing equation III. 81 with the potential used by Yakushevich in Q it is easy to see that the coupling constant 
values induced on the Yakushevich model are given by 

g = 2gt + A{l + afg,^30, K ^ 2(1 + f3f gp ^ 5.8 (III.6) 

so that the normalized Yakushevich coupling constant (the torsional one) corresponding to our choice of parameters 
is g/K ~ 5.1. 

Note that, when using the Morse expression for pairing, the constant in front of the pairing is simply the bond 
energy D = O.lkJ/mol (see eq. (III. 31) ). We decide nevertheless to divide by the same normalizing factor Ip = 
Kp{dbs + 2r)^/2 = D^^ so that all dimensionless coupling constants are unchanged except for the pairing one which 
becomes g' = 1//Lt^ ~ 0.013. 



C. Dispersion relations 

The linearized version of the composite Lagrangian is 

^Yiwi^^n,!-! ^n,ii ^n^ii ^n,i) — 



N 2 . 



n— 1 i—1 
N-1 2 



^n,t + (1 + a) 



n— 1 i—1 
N-1 2 ^ 

Yl 2^^s{dbs + r)2 [Av3„,, + (1 + a)Aen,^f 

n—1 i—1 

^ 1 

l^Kpidbs + 2r)2 [(^„ 1 + <^„ 2 + (1 + P){OnA + en,2)Y 



2 

n=l 
N 2 



n—1 i—1 



leading to the following equations of motion: 

' Iten,^ + (1 + a)Is[^n,^ + (1 + = 



Is[Vn,i + (1 + Q)'OnA 



-K,{dbs + r)2(l + a)[U^na + (1 + a)Ue,,A 

-Kp{dbs + 2r)2(l + [3) [^„,i + <^„,2 + (1 + f5){enA + 0«.2)] 

-Kb [—On+bA+1 + 26n,i — 6*71-5,1+1] 

-Ks{dbs + r)2p(p„^, + (1 + a)ne^A 

-Kp{dbs + 2rf [<^„,i + ipna + (1 + P){On,l + 6I„,2)] 



(III.7) 



(III.8 



where D^n.i = ~On+i,i + 20„.i — 9n-i,i and similarly for ip. What we are interested in are the conditions for the 
existence of wave solutions for the equations above, i.e. in the form 



q . — , . pi(kSn+ujt) . — <T> , . pi(k5n+u}t) 



J = 1,2 



In the approximation a ~ /3, namely considering the bases as pointlike, it is easy to simplify a pair of equations by 
multiplying the second equation by (1 + a) and substracting it from the first one so that we are left with the pair 



-Ktne^,, - Kh [-9n+5,,+i + 29. 



n,i "^n — b, 't+lj 1 



i] , * = 1,2 
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Summing and subtracting and imposing the wave form for the angles on the two equations above we find, by imposing 
the condition for the existence of non-trivial waves, the first pair of dispersion relations: 

ojf = 4{Kt/It)sm'^{kS/2) + 2{Kh/It)[l + cos{5kS)] 

ujI = ^{KtlIt)sa?{Ul2) ^ A,{jihlh)'&v[?{^kbl'l) (IIL9) 

The second pair of the dispersion relations can then be easily extracted from the remaining pair of equations after the 
change of coordinates ?/'n,i = </'n,i + (1 + a)^n,i, the approximation d^s + 2r ~ dhs + t (which in turn implies a — ff) 
and finally the fact that Ig/ {dbs + the mass to;, of the base; the pair of equations then become 

I mb'ipn,! = -Ksni>n,l ~ Kp [V^n^i + l/'n,2] 

Restricting the equations on the wave solutions and summing and subtracting them we find the other pair of dispersion 
relations 

ujj = 4:{K,/mb)sm'^{kS/2) + 2Kp/mb (III.IO) 

Physically, the four dispersion relations correspond to the four oscillation modes of the system in the linear regime. 
The relations involving lui and lu2 are associated with torsional oscillations of the backbone. In case of uji there is a 
threshold for the generation of the excitation originating in the helicoidal interaction, whereas the second torsional 
mode UJ2 has no threshold and is thus also of acoustical type. The relation involving uj^ describes relative oscillations 
of the two bases in the chain with respect to the neighboring bases. As wsik) — > for A: — > there is no threshold for 
the generation of these phonon mode excitations. The dispersion relation involving lu^ describes relative oscillations 
of two bases in a pair. The threshold for the generation of the excitation is now determined by the pairing interaction. 

The dispersion relations (jIII.9IIII.10[) for values of the physical parameters given in the tables U and IIIII are plotted 
in Fig. IIII Cl there we plot uj/{27rc), where c is the speed of light (we use the, in the literature widespread, convention 
of measuring frequencies in 27rc units) versus kS/2. 

The four dispersion relations take a simple form if we consider excitations with wavelength A much bigger then the 
intrapair distance, i.e A 3> (5; this corresponds to the 6-^0 limit. We have then 

^ ql, (III. 11) 

where Ca and qa (a = 1 . . .4) are, respectively, the velocity of propagation (in the limit k ^ qa) and the excitation 
threshold. They are given by 



Cl - 5 ^{Kt - 25KhWt , 91 = '^Vk^u 

C2 = 5 ^{Kt + 2 bKh)/It, 92 = 0; ^^j^^^^) 
C3 = b^jKsjmb, qz = 0; 



C4 = 5^Ks/mb, qi = yj2Kp/mb] 

Using the values of the parameters given in the tables U and IIIII we have 

Cl ~ Okm/s, qi ~ 4.5 cm~^ ; 

C2^1km/s, 92 = 0; (III 13) 

C3 ~ 3km/s, 93 = ; 

C4 ~ 3km/s, 94 ~ 32 cm^^ , 

where ci ~ comes from the fact that we are taking Kt ~ 25 Kh (see table IIIll - this of course just means that ci is at 
least an order of magnitude smaller than the other Ci , and therefore negligible) . Speeds can be converted to base per 
seconds by dividing each Ci hy 5 — 3.4A; excitation thresholds can be converted in inverse of seconds by multiplying 
each 9i by 27rc, where c is the speed of light. 

IV. NUMERICAL ANALYSIS 

As we pointed out in sec. HH the profiles of the solitons in our model are extremals of the function of 4A'' variables 
L{0n,i, ^n.i), whcrc N is the number of nodes of the chain. The interval of variations of the Lagrangian parameters 
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FIG. 5: Graph of the dispersion relations (IIII.9IIII.10P in the first Brillouin zone. We plot tJc«/(27rc) (c is the speed of light) as 
a function of k5/2. The uj\,UL>2,'^z,(jJi are represented respectively by the thick continuous, thin dashed, thin continuous and 
thick dashed line. Units are cm^^ in the vertical axis and radiants in the horizontal axis. 
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E (MJ/mol) 
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4.6 


43 
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(0,1) 


21 


4.6 


43 


29 


0.7 


144 


(1,1) 




14 


14 




1.1 


196 



TABLE IV: Energies and diameters for soliton profiles shown in fig. |6] diameters are defined as the number of consecutive 
nodes between the angles Omin = 1/10 and Omax = 2n — 9min- 

under our study is such that solitons can have a radius of more than a hundred nodes so a reasonable value of N is 
of the order of 10'^; following Yakushevich Q we mostly used the value A'' = 2000. 

Clearly the extremals - which by the way in our case turn out to be always maxima - of such a complex function 
of so many variables can only be obtained numerically; to accomplish that, again following Yakushevich, we choose 
to use the so-called conjugate gradients method. Ready-to-use implementations of this algorithm are available in the 
main numerical libraries, in particular in Numerical Recipes (NR, [27|) and in the GNU Scientific Library (GSL, 28]). 
As a double check, most of the results presented here have been produced using both implementations; all profiles 
shown in this paper were produced with GSL. 

A. Solitons profiles in the Yakushevich model 

As a generalization of the sine-Gordon solitons, also the Yakushevich solitons are "relativistic" and in particular 
admit a limit speed beyond which no solution exists. We verified numerically that, for a fixed choice of the coupling 
constants, solitons profiles corresponding to all possible speeds vary by less than 10% from each other, so that it is 
enough to show the profiles for a single value of the speed. All profiles presented here are relative to static solitons - 
in other words non-trivial equilibria position of the double chain; in this case the kynetic term is zero and the only 
parameter left in the renormalized Lagrangian is the coupling constant g of the torsion potential. 

As starting point for the extremizing algorithm, following Yakushevich, we used 

6»„ ^ = q,7r{l + tanh [M{2n - N)]} (IV.l) 

where {qi, (72) is the topological type of the soliton and M a parameter that dictates the steepness of the profile which 
we use to probe different regions of the phase space. 



A first important feature of the model is the disappearance of the discrete solitons when the coupling constant g is 
smaller than some value gmin — 10 (see table [V]), way before the natural threshold represented by when the solitons 
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FIG. 6: Profiles of the solitons with the lowest non-trivial topological numbers for the Yakushevich model in the harmonic 
approximation for pairing (left column) and with a Morse pairing potential (right). In each picture the thin line correspond 
to the physical value for the coupling constant g = 21 and the thick one, to show how the soliton profile changes with g, to 
g = 150; the continous line is relative to the angles dn,i, the dashed one to the angles 6n,2- In the (1, 1) pictures the dashed 
lines are not visible because the angles profiles for the two angles are identical. 



become so steep that jump from to 27r in a space shorter than the chain step 6. This behaviour is due to the 
impossibility of balancing between torsion and pairing when g is too small: indeed both contributions "telescope" but 
the torsion ones arc a priori bounded. A simple way to sec how this happen is to give a look to the simplest case, 
namely the one corresponding to topological numbers (1, 1), when the discrete Lagrange equations can be reduced to 
a sort of "discrete sine- Gordon" equations, namely 

gisinAOi - sin A6ii_i) = KnOi 

Summing term by term the first n equations we get 

5(sin A6l„ - sinA6lo) = X(6'„+i - On - + 9o) 

which, considering that 6o = and that for n of the order of N/2 = 1000 also 0i ~ 0, reduces to sin A0„ = {K/g)A9n- 
Clearly when g decreases the right term increases and when it becomes bigger than 1 there can be no solution anymore. 

This feature is particularly relevant for our study because the value g = 5.4 induced by the values of the coupling 
constants we extracted from literature (see previous section) unfortunately falls in the range where no soliton arise - 
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(1,0) 


(0,1) 


(1,1) 


go 


7.05 


7.05 


14.7 



TABLE V: The transition values go for solitons instability (solitons arise only for g > go) of the (p, q) solitons. 




FIG. 7: Region of the {gt,gs) plane where no solitons arise for the (0, 1) topological type. 

at least for the basic topological numbers (1, 0), (0, 1) and (1, 1) - in the harmonic approximation for pairing. Since 
nevertheless 5.4 is rather close to the boundary values for the existence of solutions and the estimate of coupling 
constants is only qualitative, for the harmonic approximation we increased the value by a factor 4. Note that though 
this behaviour disappears when a more physical potential, i.e. a Morse one, is used for pairing, leading to the 
conclusion that the feature above belongs to the harmonic version of the model but does not play any role in realistic 
models of DNA. 

In fig. [5] we show the results for the solitons corresponding to the basic topological numbers (1, 0), (0, 1) and (1, 1) 
both in the harmonic approximation for pairing (used in Yakushevich in Q) and with a Morse potential for two values 
of the coupling constant g, one corresponding to the physical value and one corresponding to a value about an order 
of magnitude above in order to enhance the detail in the solitons profile (whose diameter increases monotonically with 
g). As expected, the profiles corresponding to the Morse potential have a bigger radius, actually by about an order 
of magnitude, than the corresponding harmonic ones. In all cases the numerical solutions appear to be very robust 
with respect to the initial configuration, so that we can affirm that in all of them there is a unique extremal. 

B. Solitons profiles in the composite model 



As a generalization of the sine-Gordon model, also in this composite model there is very little difference between 
profiles corresponding to different speeds [sH : in particular, we will consider even in this case only stationary profiles. 
As starting point for the extremizing algorithm we used the natural generalization of (jIV.l[) : 

0n,^ = g.TT {1 + tauh [M {2n - N)]} , 0„^, = . 
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21 


54 
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21 


54 


1.6 


0.7 


136 


(1,1) 




63 


30 




1.1 


154 



TABLE VL Energies and diameters for soliton profiles shown in fig. [8] and [Oj diameters are defined as the number of consecutive 
nodes between the angles 9min = 1/10 and 9max = 2tt — 6min- 
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FIG. 8; Profiles of solitons corresponding to tlie smallest non-trivial topological numbers for the composite model in the 
harmonic approximation for pairing. In the (1,0) case we compare the profile with the corresponding Yakushevich one to 
evidence how close the topological component in the composite model are to their Yakushevich counterpart. In the (0, 1) case 
we compare the profiles to the ones we get by increasing the torsional part to values corresponding to g = 150 in Yakushevich 
model, once by putting the whole contribute in Vt (thin line) and then putting it wholly in Vs (thinner line) . In both cases the 
profiles diameters get bigger, as expected, but no new qualitative feature appears. In the (1, 1) case we compare the profiles 
to the ones we get by not making the contact approximation (i.e. keeping d^q 7^ 0, thinner line) or by enlarging the helicoidal 
coupling constant (thin line). Even in this case no new feature appears. 



Since the angles 4)n,i are bound to only half of a circle, the corresponding field in the continous approximation cannot 
describe a topologically non-trivial path and therefore we call these coordinates "non-topological" and there is no 
integer number associated to them, while the topological numbers qi of the 9 angles correspond exactly to the numbers 
in the simpler model. 

A further similarity between the two systems is the disappearance of the discrete soliton for too little values of the 
coupling constants, in this case gt and gs (see fig. [7] - the helicoidal term is at least an order of magnitude smaller 
than them and it can be disregarded). This is why in all profiles relative to the composite model in the harmonic 
approximation (fig. [5]) we magnified g^ by a factor 4, bringing it to g^ = 6.5; when a Morse-like potential is used 
instead (fig. (HI) solitons profiles survive even when both torsion and stacking are turned off, so in that case we keep 
using the correct value gs = 1.8. 

In this section we show the profiles of the solitons for the base topological types (1, 0), (0, 1) and (1, 1) for both the 
harmonic pairing approximation (fig. [5]) and the more physical Morse potential (fig. [5]) comparing them with some of 



14 



the several deformations that can be tested to check their stabihty. 

1. Harmonic approximation 

In fig. m first we compare the plot of the (1,0) soliton with the corresponding soliton in the Yakushevich model; 
the two graphs are within 10% from each other, testifying that the increase in complexity of the model gives us more 
features without modifying qualitatively the successful results of the model that generalizes. 

Then we compare the profiles of the (0,1) soliton with those corresponding to a torsional coupling an order of 
magnitude bigger. Since in this model there are two distinct sources for torsion, we choose to show the profiles 
corresponding to the extremal cases, namely when all torsion is concentrated in the backbone and when it is instead 
concentrated in the bases; in this case there are bigger differences between the profiles but we can say that, even after 
magnifying by an order of magnitude the torsion with respect to the pairing, the picture stays qualitatively the same. 

Finally we compare the profiles of the (1,1) soliton with those obtained by not doing the contact approximation 
(i.e. by considering the equilibrium position of bases being at distance of SA from each other) and those obtained by 
magnifying the helicoidal contribution by a two orders of magnitude. Even in this case the differences in the profiles 
turned out negligible. 

Summarizing the results, the harmonic approximation turned out to be rather solid even in the composite model 
but it keeps suffering of the same problem already detected in the Yakushevich model, namely there is no non-trivial 
discrete soliton if the torsional part is "too small" and unfortunately the physical values of the coupling constants 
seem to follow inside this non-existence zone, or at the very least very close to its boundary. 

A second concern is that the profiles of the non-topological angles, even the ones generated with higher torsional 
couplings, in correspondance to the soliton rise from to 2-k jump from about —tt/2 to 7r/2 within a single node. 
This on one side casts some doubt on the existence of a continous conterpart, leading to possible problems for their 
movement evolution on the chain, and on the other side testifies of a quite violent struggle about the points ±7r/2, 
where two very strong non-linear forces (the pairing and the sugar wall) compete with each other; this behaviour is 
quite unwelcome first because it is unphysical, since after the ionic hydrogen bonds break (and this happens when 
they get a few A apart) they do not interact anymore, and second because such a violent non-linear interaction most 
likely destabilizes the system. 



2. Morse potential 

In fig. m first we compare the plot of the (1,0) soliton with the corresponding soliton in the Yakushevich model; 
even with the Morse potential the similarity of the profiles obtained for the two models is striking, keeping within 
10% from each other. 

Then we compare the plot of the (0, 1) soliton with the ones obtained by decreasing the well width parameter i.e. 
enlarging the well from lA to 4A. As expected we get a norrower profile since in the limit /i ^ the Morse potential 
reduces to its harmonic approximation. 

Finally we show the (1,1) soliton profiles and compare them with the ones obtained by neglecting altogether both 
torsional couplings. The profile is of course narrower but it is still non-trivial, as opposite to the case of the harmonic 
approximation when even at the physical values for the coupling constants the only profiles we obtain are the constant 
ones. 

Summarizing the numerical analysis relative to the Morse potential, we get the same nice features found in the 
harmonic approximation but we loose its worst defects. It seems rather clear therefore that a serious investigation of 
DNA's rotational dynamical properties cannot avoid using a Morse-like potential for pairing, i.e. the pairing coupling 
must be turned completely off after the distance between bases increases by a few A. 

V. CONCLUSIONS 

The numerical analysis we have performed shows the existence of solitonic solutions of our composite DNA model. 
The profiles of the topological solitons - in particular, the part relating to the topological degree of freedom - of our 
model are both qualitatively and quantitatively very similar to those of the Yakushevich model. This means that 
the most relevant (for DNA transcription) and characterizing feature of the nonlinear DNA dynamics present in the 
Yakushevich model is preserved by considering geometrically more complex and hence more realistic DNA models. 
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FIG. 9: Profiles of solitons corresponding to the smallest non-trivial topological numbers for the composite model with a Morse 

potential for pairing. In the (1,0) case wc compare the profile with the corresponding Yakuslievicli one to evidence how close 
are even in this case the topolgical component in the composite model to their Yakuslievicli counterpart. In the (0, 1) case wc 
compare the profiles to the ones we get by decreasing the well width normalized parameter ^{dbs + 2r) from 8.8 to 2 (thinner 
line). As expected, the profile with a smaller n gets narrower, due to the fact that the radius of the Morse hole got bigger and 
therefore the pairing interaction looks more like a globally harmonic one. In the (1, 1) case we compare the profiles to the ones 
we get by disregarding completely both torsional couplings. Differently from what happens in the harmonic approximation, 
when solitons become trivial for small torsions, the profile has the same qualitative features of the one corresponding to the 
physical coupling constants. 



Moreover, the topological soliton profiles of our model seem to change very little when either the physical parameters 
change in a reasonable range or also the form of the potential modelling the pairing interaction is modified to a 
more realistic form with a sole exception, namely the replacement of the pairing Morse potential with its harmonic 
approximation, which works fine close to the equilibrium position but fails badly when the bases flip. 

In particular, the forms of the topological solitons are very little sensitive to the interchange of torsional and 
stacking coupling constant. This feature adds other reasons why the Yakushevich model, although based on a strong 
simplification of the DNA geometry, works quite well in describing solitonic excitations. The Yakushevich model, 
indeed, does not distinguish between torsional and stacking interaction; but, as we have shown, this distinction is 
not relevant - at least as long as one is only interested in the existence and form of the soliton solutions. The 
"compositeness" of our model becomes relevant - and rather crucial - when it comes on the one hand to allowing the 
existence of solitons together with requiring a physically realistic choice of the physical parameters characterizing the 
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DNA, and on the other hand to have also predictions fitting experimental observations for what concerns quantities 
related to small amplitude dynamics, such as transverse phonons speed. In other words, the somewhat more detailed 
description of DNA dynamics provided by our model allows it to be effective - with the same parameters - across 
regimes, and provide meaningful quantities in both the linear and the fully nonlinear regime. 

We expect that the model considered here is the simplest DNA model describing rotational degrees of freedom 
which, with physically realistic values of the coupling constants and other parameters, allows for the existence of 
topological solitons and at the same time is also compatible with observed values of bound energies and phonon 
speeds in DNA. We also expect that solitons may move for considerable distances in this model even in presence of 
realistic inhomogeneities thanks to the fact that the component that supports the solitons motion, i.e. the sugar- 
phosphate group, is homogeneous and therefore by separating it from the bases we can consider the inhomogeneities 
as a perturbative effect. We leave to a future paper the study of inhomogeneities and time evolution in this model. 
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